Notebook Imports¶
In [1]:
from sklearn.datasets import load_boston
from sklearn.model_selection import train_test_split
from sklearn.linear_model import LinearRegression
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns
import numpy as np
import statsmodels.api as sm
from statsmodels.stats.outliers_influence import variance_inflation_factor
%matplotlib inline
Gather Data¶
In [2]:
boston_dataset = load_boston()
In [3]:
type(boston_dataset)
Out[3]:
In [4]:
boston_dataset
Out[4]:
In [5]:
dir(boston_dataset)
Out[5]:
In [6]:
print(boston_dataset.DESCR)
Data points and features¶
In [7]:
type(boston_dataset.data)
Out[7]:
In [8]:
boston_dataset.data.shape # chaining dot notation
Out[8]:
In [9]:
boston_dataset.feature_names
Out[9]:
In [10]:
# Actual prices in thousands (000s)
boston_dataset.target
Out[10]:
Data exploration with Pandas dataframes¶
In [11]:
# Create a pandas dataframe
data = pd.DataFrame(data=boston_dataset.data, columns=boston_dataset.feature_names)
# Add column with the price (target)
data['PRICE'] = boston_dataset.target
In [12]:
data.head() # The top rows look like this
Out[12]:
In [13]:
data.tail() # Rows at bottom of dataframe look like this:
Out[13]:
In [14]:
data.count() # show us the number of row
Out[14]:
Cleaning data - check for missing values¶
In [15]:
pd.isnull(data).any()
Out[15]:
In [16]:
data.info()
Visualising Data - Histograms, Distributions and Bar Charts¶
In [17]:
plt.figure(figsize=(10, 6))
plt.hist(data['PRICE'], bins=50, ec='black', color='#2196f3')
plt.xlabel('Price in 000s')
plt.ylabel('Nr. of Houses')
plt.show()
In [18]:
plt.figure(figsize=(10, 6))
sns.distplot(data['PRICE'], bins=50, hist=True, kde=False, color='#fbc02d')
plt.show()
In [19]:
plt.figure(figsize=(10, 6))
plt.hist(data['RM'], ec='black', color='#00796b')
plt.xlabel('Average Number of Rooms')
plt.ylabel('Nr. of Houses')
plt.show()
In [20]:
data['RM'].mean()
Out[20]:
In [21]:
# Challenge: Create a meaningful histogram for RAD using matplotlib ... in royal purple
plt.figure(figsize=(10, 6))
plt.hist(data['RAD'], bins=24, ec='black', color='#7b1fa2', rwidth=0.5)
plt.xlabel('Accessibility to Highways')
plt.ylabel('Nr. of Houses')
plt.show()
In [22]:
data['RAD'].value_counts()
Out[22]:
In [23]:
frequency = data['RAD'].value_counts()
#type(frequency)
#frequency.index
#frequency.axes[0]
plt.figure(figsize=(10, 6))
plt.xlabel('Accessibility to Highways')
plt.ylabel('Nr. of Houses')
plt.bar(frequency.index, height=frequency)
plt.show()
In [24]:
data['CHAS'].value_counts()
Out[24]:
Descriptive Statistics¶
In [25]:
data['PRICE'].min()
Out[25]:
In [26]:
data['PRICE'].max()
Out[26]:
In [27]:
data.min()
Out[27]:
In [28]:
data.max()
Out[28]:
In [29]:
data.mean()
Out[29]:
In [30]:
data.median()
Out[30]:
In [31]:
data.describe()
Out[31]:
In [32]:
data['PRICE'].corr(data['RM'])
Out[32]:
In [33]:
# Challenge: Calculate the correlation between property prices and the pupil teacher ratio
data['PRICE'].corr(data['PTRATIO'])
Out[33]:
In [34]:
data.corr() # Pearson Correlation Coefficients
Out[34]:
In [35]:
mask = np.zeros_like(data.corr())
triangle_indices = np.triu_indices_from(mask)
mask[triangle_indices] = True
mask
Out[35]:
In [36]:
plt.figure(figsize=(16,10))
sns.heatmap(data.corr(), mask=mask, annot=True, annot_kws={"size": 14})
sns.set_style('white')
plt.xticks(fontsize=14)
plt.yticks(fontsize=14)
plt.show()
In [37]:
# Challenge: Picture the relationship between pollution and distance in your head
# Then create a scatter plot between DIS and NOX
In [38]:
nox_dis_corr = round(data['NOX'].corr(data['DIS']), 3)
plt.figure(figsize=(9, 6), dpi=300)
plt.scatter(x=data['DIS'], y=data['NOX'], alpha=0.6, s=80, color='indigo')
plt.title(f'DIS vs NOX (Correlation {nox_dis_corr})', fontsize=14)
plt.xlabel('DIS - Distance from employment', fontsize=14)
plt.ylabel('NOX - Nitric Oxide Pollution', fontsize=14)
plt.show()
In [ ]:
sns.set()
sns.set_context('talk')
sns.set_style('whitegrid')
sns.jointplot(x=data['DIS'], y=data['NOX'], height=7, color='indigo', joint_kws={'alpha':0.5})
plt.show()
In [40]:
sns.set()
sns.set_context('talk')
sns.set_style('whitegrid')
sns.jointplot(x=data['DIS'], y=data['NOX'], kind='hex', size=7, color='blue')
plt.show()
In [41]:
sns.set()
sns.set_context('talk')
sns.set_style('whitegrid')
sns.jointplot(x=data['TAX'], y=data['RAD'], size=7, color='darkred', joint_kws={'alpha':0.5})
plt.show()
In [ ]:
sns.lmplot(x='TAX', y='RAD', data=data, height=7)
plt.show()
In [43]:
# Challenge: Create a scatter plot between the house prices and the number of rooms (RM)
In [44]:
rm_tgt_corr = round(data['RM'].corr(data['PRICE']), 3)
plt.figure(figsize=(9, 6))
plt.scatter(x=data['RM'], y=data['PRICE'], alpha=0.6, s=80, color='skyblue')
plt.title(f'RM vs PRICE (Correlation {rm_tgt_corr})', fontsize=14)
plt.xlabel('RM - Median nr of rooms', fontsize=14)
plt.ylabel('PRICE - property price in 000s', fontsize=14)
plt.show()
In [45]:
sns.lmplot(x='RM', y='PRICE', data=data, size=7)
plt.show()
In [46]:
%%time
sns.pairplot(data)
plt.show()
In [47]:
%%time
sns.pairplot(data, kind='reg', plot_kws={'line_kws':{'color': 'cyan'}})
plt.show()
Training & Test Dataset Split¶
In [87]:
prices = data['PRICE']
features = data.drop('PRICE', axis=1)
X_train, X_test, y_train, y_test = train_test_split(features, prices,
test_size=0.2, random_state=10)
# % of training set
len(X_train)/len(features)
Out[87]:
In [88]:
# % of test data set
X_test.shape[0]/features.shape[0]
Out[88]:
Multivariable Regression¶
In [89]:
regr = LinearRegression()
regr.fit(X_train, y_train)
#Challenge: print out r-squared for training and test datasets
print('Training data r-squared:', regr.score(X_train, y_train))
print('Test data r-squared:', regr.score(X_test, y_test))
print('Intercept', regr.intercept_)
pd.DataFrame(data=regr.coef_, index=X_train.columns, columns=['coef'])
Out[89]:
Data Transformations¶
In [90]:
data['PRICE'].skew()
Out[90]:
In [99]:
y_log = np.log(data['PRICE'])
y_log.tail()
Out[99]:
In [100]:
y_log.skew()
Out[100]:
In [93]:
sns.distplot(y_log)
plt.title(f'Log price with skew {y_log.skew()}')
plt.show()
In [94]:
sns.lmplot(x='LSTAT', y='PRICE', data=data, size=7,
scatter_kws={'alpha': 0.6}, line_kws={'color':'darkred'})
plt.show()
In [96]:
transformed_data = features
transformed_data['LOG_PRICE'] = y_log
sns.lmplot(x='LSTAT', y='LOG_PRICE', data=transformed_data, height=7,
scatter_kws={'alpha': 0.6}, line_kws={'color':'cyan'})
plt.show()
Regression using log prices¶
In [97]:
prices = np.log(data['PRICE']) # Use log prices
features = data.drop('PRICE', axis=1)
X_train, X_test, y_train, y_test = train_test_split(features, prices,
test_size=0.2, random_state=10)
regr = LinearRegression()
regr.fit(X_train, y_train)
print('Training data r-squared:', regr.score(X_train, y_train))
print('Test data r-squared:', regr.score(X_test, y_test))
print('Intercept', regr.intercept_)
pd.DataFrame(data=regr.coef_, index=X_train.columns, columns=['coef'])
Out[97]:
In [58]:
# Charles River Property Premium
np.e**0.080475
Out[58]:
p values & Evaluating Coefficients¶
In [59]:
X_incl_const = sm.add_constant(X_train)
model = sm.OLS(y_train, X_incl_const)
results = model.fit()
#results.params
#results.pvalues
pd.DataFrame({'coef': results.params, 'p-value': round(results.pvalues, 3)})
Out[59]:
Testing for Multicollinearity¶
$$ TAX = \alpha _0 + \alpha _1 RM + \alpha _2 NOX + ... + \alpha _{12}LSTAT $$$$ VIF _{TAX} = \frac{1}{(1 - R _{TAX} ^ 2)} $$In [60]:
variance_inflation_factor(exog=X_incl_const.values, exog_idx=1)
#type(X_incl_const)
Out[60]:
In [61]:
# Challenge: print out the number of columns in X_incl_const
len(X_incl_const.columns)
X_incl_const.shape[1]
Out[61]:
In [62]:
# Challenge: write a for loop that prints out all the VIFs for all the features
for i in range(X_incl_const.shape[1]):
print(variance_inflation_factor(exog=X_incl_const.values, exog_idx=i))
print('All done!')
In [63]:
vif = [] # empty list
for i in range(X_incl_const.shape[1]):
vif.append(variance_inflation_factor(exog=X_incl_const.values, exog_idx=i))
print(vif)
In [64]:
vif = [variance_inflation_factor(exog=X_incl_const.values,
exog_idx=i) for i in range(X_incl_const.shape[1])]
pd.DataFrame({'coef_name': X_incl_const.columns,
'vif': np.around(vif, 2)})
Out[64]:
In [65]:
import this
Model Simplification & the BIC¶
In [66]:
# Original model with log prices and all features
X_incl_const = sm.add_constant(X_train)
model = sm.OLS(y_train, X_incl_const)
results = model.fit()
org_coef = pd.DataFrame({'coef': results.params, 'p-value': round(results.pvalues, 3)})
# Challenge: find and check official docs for results object and print out BIC & r-squared
print('BIC is', results.bic)
print('r-squared is', results.rsquared)
In [67]:
# Reduced model #1 excluding INDUS
X_incl_const = sm.add_constant(X_train)
X_incl_const = X_incl_const.drop(['INDUS'], axis=1)
model = sm.OLS(y_train, X_incl_const)
results = model.fit()
coef_minus_indus = pd.DataFrame({'coef': results.params, 'p-value': round(results.pvalues, 3)})
print('BIC is', results.bic)
print('r-squared is', results.rsquared)
In [68]:
# Reduced model #2 excluding INDUS and AGE
X_incl_const = sm.add_constant(X_train)
X_incl_const = X_incl_const.drop(['INDUS', 'AGE'], axis=1)
model = sm.OLS(y_train, X_incl_const)
results = model.fit()
reduced_coef = pd.DataFrame({'coef': results.params, 'p-value': round(results.pvalues, 3)})
print('BIC is', results.bic)
print('r-squared is', results.rsquared)
In [69]:
frames = [org_coef, coef_minus_indus, reduced_coef]
pd.concat(frames, axis=1)
Out[69]:
Residuals & Residual Plots¶
In [70]:
# Modified model: transformed (using log prices) & simplified (dropping two features)
prices = np.log(data['PRICE']) # Use log prices
features = data.drop(['PRICE', 'INDUS', 'AGE'], axis=1)
X_train, X_test, y_train, y_test = train_test_split(features, prices,
test_size=0.2, random_state=10)
# Using Statsmodel
X_incl_const = sm.add_constant(X_train)
model = sm.OLS(y_train, X_incl_const)
results = model.fit()
# Residuals
# residuals = y_train - results.fittedvalues
# results.resid
# Graph of Actual vs. Predicted Prices
corr = round(y_train.corr(results.fittedvalues), 2)
plt.scatter(x=np.e**y_train, y=np.e**results.fittedvalues, c='blue', alpha=0.6)
plt.plot(np.e**y_train, np.e**y_train, color='cyan')
plt.xlabel('Actual prices 000s $y _i$', fontsize=14)
plt.ylabel('Prediced prices 000s $\hat y _i$', fontsize=14)
plt.title(f'Actual vs Predicted prices: $y _i$ vs $\hat y_i$ (Corr {corr})', fontsize=17)
plt.show()
# Residuals vs Predicted values
plt.scatter(x=results.fittedvalues, y=results.resid, c='navy', alpha=0.6)
plt.xlabel('Predicted log prices $\hat y _i$', fontsize=14)
plt.ylabel('Residuals', fontsize=14)
plt.title('Residuals vs Fitted Values', fontsize=17)
plt.show()
# Mean Squared Error & R-Squared
reduced_log_mse = round(results.mse_resid, 3)
reduced_log_rsquared = round(results.rsquared, 3)
In [71]:
# Distribution of Residuals (log prices) - checking for normality
resid_mean = round(results.resid.mean(), 3)
resid_skew = round(results.resid.skew(), 3)
sns.distplot(results.resid, color='navy')
plt.title(f'Log price model: residuals Skew ({resid_skew}) Mean ({resid_mean})')
plt.show()
In [72]:
# Challenge: Using the original model with all the features and normal prices generate:
# Plot of actual vs predicted prices (incl. correlation) using a different colour
# Plot of residuals vs. predicted prices
# Plot of distribution of residuals (incl. skew)
# Analyse the results.
In [73]:
# Original model: normal prices & all features
prices = data['PRICE']
features = data.drop(['PRICE'], axis=1)
X_train, X_test, y_train, y_test = train_test_split(features, prices,
test_size=0.2, random_state=10)
X_incl_const = sm.add_constant(X_train)
model = sm.OLS(y_train, X_incl_const)
results = model.fit()
# Graph of Actual vs. Predicted Prices
corr = round(y_train.corr(results.fittedvalues), 2)
plt.scatter(x=y_train, y=results.fittedvalues, c='indigo', alpha=0.6)
plt.plot(y_train, y_train, color='cyan')
plt.xlabel('Actual prices 000s $y _i$', fontsize=14)
plt.ylabel('Prediced prices 000s $\hat y _i$', fontsize=14)
plt.title(f'Actual vs Predicted prices: $y _i$ vs $\hat y_i$ (Corr {corr})', fontsize=17)
plt.show()
# Residuals vs Predicted values
plt.scatter(x=results.fittedvalues, y=results.resid, c='indigo', alpha=0.6)
plt.xlabel('Predicted prices $\hat y _i$', fontsize=14)
plt.ylabel('Residuals', fontsize=14)
plt.title('Residuals vs Fitted Values', fontsize=17)
plt.show()
# Residual Distribution Chart
resid_mean = round(results.resid.mean(), 3)
resid_skew = round(results.resid.skew(), 3)
sns.distplot(results.resid, color='indigo')
plt.title(f'Residuals Skew ({resid_skew}) Mean ({resid_mean})')
plt.show()
# Mean Squared Error & R-Squared
full_normal_mse = round(results.mse_resid, 3)
full_normal_rsquared = round(results.rsquared, 3)
In [74]:
# Model Omitting Key Features using log prices
prices = np.log(data['PRICE'])
features = data.drop(['PRICE', 'INDUS', 'AGE', 'LSTAT', 'RM', 'NOX', 'CRIM'], axis=1)
X_train, X_test, y_train, y_test = train_test_split(features, prices,
test_size=0.2, random_state=10)
X_incl_const = sm.add_constant(X_train)
model = sm.OLS(y_train, X_incl_const)
results = model.fit()
# Graph of Actual vs. Predicted Prices
corr = round(y_train.corr(results.fittedvalues), 2)
plt.scatter(x=y_train, y=results.fittedvalues, c='#e74c3c', alpha=0.6)
plt.plot(y_train, y_train, color='cyan')
plt.xlabel('Actual log prices $y _i$', fontsize=14)
plt.ylabel('Predicted log prices $\hat y _i$', fontsize=14)
plt.title(f'Actual vs Predicted prices with omitted variables: $y _i$ vs $\hat y_i$ (Corr {corr})', fontsize=17)
plt.show()
# Residuals vs Predicted values
plt.scatter(x=results.fittedvalues, y=results.resid, c='#e74c3c', alpha=0.6)
plt.xlabel('Predicted prices $\hat y _i$', fontsize=14)
plt.ylabel('Residuals', fontsize=14)
plt.title('Residuals vs Fitted Values', fontsize=17)
plt.show()
# Mean Squared Error & R-Squared
omitted_var_mse = round(results.mse_resid, 3)
omitted_var_rsquared = round(results.rsquared, 3)
In [75]:
pd.DataFrame({'R-Squared': [reduced_log_rsquared, full_normal_rsquared, omitted_var_rsquared],
'MSE': [reduced_log_mse, full_normal_mse, omitted_var_mse],
'RMSE': np.sqrt([reduced_log_mse, full_normal_mse, omitted_var_mse])},
index=['Reduced Log Model', 'Full Normal Price Model', 'Omitted Var Model'])
Out[75]:
In [76]:
# Challenge: Our estimate for a house price is $30,000. Calculate the upper and lower bound
# for a 95% prediction interval using the reduced log model
print('1 s.d. in log prices is', np.sqrt(reduced_log_mse))
print('2 s.d. in log prices is', 2*np.sqrt(reduced_log_mse))
upper_bound = np.log(30) + 2*np.sqrt(reduced_log_mse)
print('The upper bound in log prices for a 95% prediction interval is ', upper_bound)
print('The upper bound in normal prices is $', np.e**upper_bound * 1000)
lower_bound = np.log(30) - 2*np.sqrt(reduced_log_mse)
print('The lower bound in log prices for a 95% prediction interval is ', lower_bound)
print('The lower bound in normal prices is $', np.e**lower_bound * 1000)